/* ====================================================
#   Copyright (C) 2021-2025  All rights reserved.
#
#   Author        : Jinjing Li
#   Email         : josiahleeoaa@outlook.com
#   File Name     : libAlphaN.h
#   Last Modified : 2025-05-05
#   Describe      : Library for (alpha,n) reaction neutron generation
#                   - Updated parameter lists (removed unused params)
#                   - Added detailed API documentation
#
#   Features:
#   - Provides functions for alpha-induced neutron generation
#   - Supports multiple data versions and materials
#   - Includes gamma de-excitation calculations
#   - Uses NEUCBOT nuclear reaction code outputs
#
#   Dependencies:
#   - ROOT (for histogram handling)
#   - CLHEP (for physical constants)
#
# ====================================================*/

#ifndef _LIBALPHAN_H
#define _LIBALPHAN_H

#define MAX_LEVELS 30

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <cmath>
#include <vector>
#include <map>
#include <algorithm>
#include <iostream>
#include <CLHEP/Units/PhysicalConstants.h>

#include "TFile.h"
#include "TSystem.h"
#include "TH1D.h"
#include "TRandom3.h"

// Path to alpha-n reaction data files
static std::string alpha_data_path = "./inputs/alphaN/";

// External variables (defined in main program)
extern std::string pwd;                 // Executable directory path
extern std::string data_version_global; // Selected data version

/**
 * @brief 3D direction vector structure with spherical coordinates
 *
 * Stores both spherical (rho,theta,phi) and Cartesian (x,y,z) coordinates
 * with automatic conversion between representations.
 */
struct dirs
{
    double rho;    ///< Radial distance (magnitude)
    double theta;  ///< Polar angle (radians)
    double phi;    ///< Azimuthal angle (radians)
    double x_norm; ///< Normalized x component
    double y_norm; ///< Normalized y component
    double z_norm; ///< Normalized z component

    /**
     * @brief Construct from Cartesian velocity vector [x,y,z]
     *
     * Normalized components will be fractions in each direction
     *
     * @param arr Pointer to 3-element array of Cartesian coordinates
     */
    dirs(double *arr)
    {
        rho = sqrt(pow(arr[0], 2) + pow(arr[1], 2) + pow(arr[2], 2));
        double z = arr[2];
        double radius = sqrt(pow(arr[0], 2) + pow(arr[1], 2));
        theta = acos(z / rho);
        phi = acos(arr[0] / radius);
        x_norm = arr[0] / rho;
        y_norm = arr[1] / rho;
        z_norm = arr[2] / rho;
    }

    /**
     * @brief Construct from spherical coordinates
     *
     * Normalized components will be absolute values in each direction
     *
     * @param rho_ Radial distance
     * @param theta_ Polar angle (radians)
     * @param phi_ Azimuthal angle (radians)
     */
    dirs(double rho_, double theta_, double phi_)
    {
        rho = rho_;
        theta = theta_;
        phi = phi_;
        x_norm = rho * sin(theta) * cos(phi);
        y_norm = rho * sin(theta) * sin(phi);
        z_norm = rho * cos(theta);
    }
};

/**
 * @brief Simple object structure for gamma decay candidates
 *
 * Stores gamma energy and decay mode information
 */
struct simpleObj
{
    double energy;        ///< Gamma ray energy (MeV)
    std::string key_mode; ///< Decay mode identifier
};

// ==============================================
// Main Generator Functions (Top Level Interface)
// ==============================================

/**
 * @brief Generate neutrons from (alpha,n) reactions for single excited state
 *
 * @param isotope Target isotope number (e.g. 90232 for Th-232)
 * @param target Target nucleus (e.g. "b11" for Boron-11)
 * @param state_idx Excitation state index (0=ground state, -1=all states)
 * @param nEvents Number of neutron events to generate
 * @param stream Output file stream (stdout or file pointer)
 * @param database Nuclear database name ("JENDL", "ENDF", etc.)
 *
 * @details
 * This function generates neutron events from (alpha,n) reactions using:
 * 1. NEUCBOT nuclear reaction code outputs
 * 2. Isotropic angular distribution sampling
 * 3. Optional gamma de-excitation simulation
 *
 * Physics model:
 * - Neutron energy sampled from pre-calculated spectra
 * - Angular distribution assumed isotropic in CM frame
 * - Gamma cascades simulated for excited states
 *
 * Output format:
 * - Each event starts with particle count
 * - Followed by particle data lines (PDG code, px, py, pz, E)
 */
void generator_alphan_basedon_sundean(int isotope, const std::string &target, int state_idx,
                                      unsigned int nEvents, FILE *stream,
                                      const std::string &database);

/**
 * @brief Generate neutrons from (alpha,n) reactions for multiple states
 *
 * @param isotope Target isotope number
 * @param target Target nucleus
 * @param nEvents Number of events to generate
 * @param stream Output file stream
 * @param Material Target material
 * @param database Nuclear database name
 * @param nums Number of states to consider
 */
void generator_alphan_basedon_sundean(int isotope, const std::string &target,
                                      unsigned int nEvents, FILE *stream,
                                      const std::string &database);

// ==============================================
// Core Processing Functions
// ==============================================

/**
 * @brief Get neutron energies from NEUCBOT output
 *
 * @param key Reaction key
 * @param Material Target material
 * @param nEvents Number of events to generate
 * @param neutron_energy_vals Output vector for neutron energies
 * @param state_idx Excitation state index
 */
void get_neutron_energy_from_sundean(const std::string &key,
                                     unsigned int nEvents, std::vector<double> &neutron_energy_vals);

/**
 * @brief Get gamma ray energies for single state
 *
 * @param target Target nucleus
 * @param material Target material
 * @param isotope Isotope number
 * @param state_idx Excitation state index
 * @param key_mode Output vector for gamma modes
 * @param intensity Output vector for gamma intensities
 * @param energy Output vector for gamma energies
 */
void get_gamma_energy(const std::string &target,
                      int state_idx, std::vector<std::string> &key_mode,
                      std::vector<double> &intensity, std::vector<double> &energy);

/**
 * @brief Get gamma ray energies for multiple states
 *
 * @param target Target nucleus
 * @param Material Target material
 * @param isotope Isotope number
 * @param database Nuclear database name
 * @param key_mode Array of output vectors for gamma modes
 * @param intensity Array of output vectors for gamma intensities
 * @param energy Array of output vectors for gamma energies
 * @param branch_ratios Output vector for branch ratios
 * @param num_states Number of states to consider
 */
void get_gamma_energy(const std::string &target, int isotope,
                      const std::string &database, std::vector<std::string> key_mode[],
                      std::vector<double> intensity[], std::vector<double> energy[],
                      std::vector<double> &branch_ratios);

/**
 * @brief Generate de-excitation gamma rays
 *
 * @param key_mode Vector of gamma modes
 * @param intensity Vector of gamma intensities
 * @param energy Vector of gamma energies
 * @param Eout Output vector for generated gamma energies
 */
void generate_deexcitation_gamma(const std::vector<std::string> &key_mode,
                                 const std::vector<double> &intensity,
                                 const std::vector<double> &energy,
                                 std::vector<double> &Eout);

// ==============================================
// Utility Functions
// ==============================================

/**
 * @brief Get source filename for alpha-n data based on key
 *
 * @param key Unique identifier for the reaction
 * @param Material Target material name
 * @param data_version Data version string (optional)
 * @return std::string Full path to source data file
 */
std::string get_alphan_source_filename_from_key(const std::string &key);

/**
 * @brief Get ROOT filename for alpha-n data based on key
 *
 * @param key Unique identifier for the reaction
 * @param Material Target material name
 * @param data_version Data version string (optional)
 * @return std::string Full path to ROOT data file
 */
std::string get_alphan_root_filename_from_key(const std::string &key);

/**
 * @brief Generate unique key for alpha-n reaction
 *
 * @param isotope Target isotope number
 * @param Material Target material name
 * @param target Target nucleus (e.g. "b11")
 * @param state_idx Excitation state index
 * @param database Nuclear database name
 * @return std::string Generated reaction key
 */
std::string get_alphan_key(int isotope, const std::string &target,
                           int state_idx, const std::string &database);

/**
 * @brief Get filename for target nucleus data
 *
 * @param target Target nucleus (e.g. "b11")
 * @param state_idx Excitation state index
 * @return std::string Target data filename
 */
std::string get_alphan_target_filename(const std::string &target, int state_idx);

/**
 * @brief Read NEUCBOT nuclear reaction code output
 *
 * @param input_filename Input filename from NEUCBOT
 * @param output_filename Output filename to write processed data
 * @return int Status code (0 = success)
 */
int read_sundean_output(const char *input_filename, const char *output_filename);

bool directoryExists(const std::string &path);

bool isRegularFile(const std::string &path);

bool createDirectory(const std::string &path);

void TestDirectory(const std::string &source_name);

#endif // _LIBALPHAN_H
